home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Best of www.BestZips.com (Collector's Edition)
/
Best of WWW.BESTZIPS.COM Collector's Edition (JCSM Shareware) (JCS Marketing).ISO
/
prgtools
/
tn2.zip
/
FFT.T
< prev
next >
Wrap
Text File
|
1996-11-15
|
3KB
|
157 lines
%
% "fft.t" computes Fast Fourier Transform
%
% Sample program for the T Interpreter by:
%
% Stephen R. Schmitt
% 962 Depot Road
% Boxborough, MA 01719
%
const NU : int := 7
const DIM : int := 2^NU % number of sample points
const T : real := 0.125 % sample interval
type carray : array[DIM] of complex
const PI : real := 2 * arccos( 0 )
var W : carray
program
var k, p, n : int
var f, s, t : real
var c : complex
var x : carray % time series
put "time series:"
put "sec ampl - 0 +"
for k := 0...DIM-1 do
s := sin( 4 * k * T )
t := k * T
put t:6:2, repeat( " ", round( 10 * ( 1 + s ) ) ) & "*"
x[k].re := s
x[k].im := 0
end for
put "\ncrunching fft . . .\n"
fft( x )
put "fft spectrum:"
put "freq power"
p := DIM div 2
for n := p...DIM-1 do
c.re := x[n].re
c.im := x[n].im
f := ( n - DIM ) / ( DIM * T )
put f:6:2, repeat( " ", round( cabs( c ) ) ) & "*"
end for
for n := 0...p do
c.re := x[n].re
c.im := x[n].im
f := n / ( DIM * T )
put f:6:2, repeat( " ", round( cabs( c ) ) ) & "*"
end for
end program
%
% fast fourier transform
%
procedure fft( var x : carray )
var t : complex
var a, c, s : real
var h, i, j, k, m, n2, nu1, p : int
label another :
n2 := DIM div 2
nu1 := NU - 1
k := 0
for h := 1...NU do
another:
for i := 1...n2 do
m := k div 2^nu1
p := bitrev( m )
a := 2 * PI * p / DIM
c := cos( a )
s := sin( a )
j := k + n2
t.re := x[j].re * c + x[j].im * s
t.im := x[j].im * c - x[j].re * s
x[j].re := x[k].re - t.re
x[j].im := x[k].im - t.im
x[k].re := x[k].re + t.re
x[k].im := x[k].im + t.im
incr k
end for
k := k + n2
if k < DIM then
goto another
end if
k := 0
decr nu1
n2 := n2 div 2
end for
for k := 0...DIM-1 do
i := bitrev( k )
if i > k then
t.re := x[k].re
t.im := x[k].im
x[k].re := x[i].re
x[k].im := x[i].im
x[i].re := t.re
x[i].im := t.im
end if
end for
end procedure
%
% compute bit reversed index for fft
%
function bitrev( j : int ) : int
var i, j1, j2, rv : int
j1 := j
rv := 0
for i := 1...NU do
j2 := j1 div 2
rv := rv * 2 + ( j1 - 2 * j2 )
j1 := j2
end for
return rv
end function